Orbital motion of spiral waves in excitable media 
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Spiral waves in active media react to small perturbations as particle- like objects. Here we apply the 
asymptotic theory to the interaction of spiral waves with a localized inhomogeneity, which leads to 
a novel prediction: drift of the spiral rotation centre along circular orbits around the inhomogeneity. 
The stationary orbits have alternating stability and fixed radii, determined by the properties of 
the bulk medium and the type of inhomogeneity, while the drift speed along an orbit depends on 
the strength of the inhomogeneity. Direct simulations confirm the validity and robustness of the 
theoretical predictions and show that these unexpected effects should be observable in experiment. 

PACS numbers: 05.45.-a, 82.40.Ck, 87.18. Hf, 87.19. Hh 



The interest in the dynamics of spiral waves as regimes 
of self-organization has considerably broadened in the 
last decades, as they have been found in ever more phys- 
ical systems of diverse types (magnetic films [1], liq- 
uid crystals [2], nonlinear optics [3], new chemical sys- 
tems [4], and in population [5], tissue [6], and subcel- 
lular [7] biology). In a perfectly uniform medium the 
core of a spiral wave may be anywhere, depending on 
initial conditions. However, real systems are always het- 
erogeneous, and therefore spiral drift due to inhomogene- 
ity is of great practical interest to applications. Under- 
standably, such drift has been mostly studied in excitable 
chemical reactions and the heart, where drift due to a 
gradient of medium properties [8, 9] and pinning [29] (an- 
choring, trapping) to a localized inhomogeneity [10-12] 
have been observed in experiments and simulations. In- 
teraction with localized inhomogeneity can be considered 
to be a particular case of the general phenomenon of vor- 
tex pinning to material defects [13]. 

Here we identify a new type of spiral wave dynam- 
ics: precession around a localized inhomogeneity along a 
stable circular orbit. We predict this novel phenomenon 
theoretically, describe its key features, and confirm it by 
numerical simulations. We argue that this orbital move- 
ment of spiral waves is robust and prevalent, has non- 
trivial and surprising consequences for applications and 
should be directly observable in experiments. 

We consider reaction-diffusion equations, which is the 
most popular class of models describing spiral waves: 



c)tu = f (u, p) + DV^u, 



(1) 



where u,f G M% D e W""', p e K", u(f,t) is the dy- 
namic vector field, f £ M^, p(r) = po + Pi{r), jpij ^ 1, 
is the vector of parameters, D is diffusion matrix. For 
P = Po = const, system (1) is assumed to have spiral 



wave solutions rotating with angular velocity cj (taken 
here to be clockwise for a; > 0), 



u = U(/9,i9 + u;i-$), 



(2) 



where (p, i?) arc polar coordinates defined with respect to 
the center of rotation R = (X, y)""", and $ is the initial 
rotation phase. 

In the presense of a small perturbation pi(7^) 7^ 0, the 
spiral's center of rotation R = X + iY is not constant but 
slowly evolves with the equation of motion 
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[W (p,6')]+h(r,T) dVdr, 



t — Tv/u. 



where p = p{r — R) and 6 = 'd{f— R) + cur - 



(3) 

$ are polar 



coordinates in the corotating frame of reference, and h is 
the perturbation to the right-hand side of Eq. (1). Func- 
tion W is called the response function (RF) and defines 
the sensitivity of the spiral wave position with respect to 
perturbations in different places. Technically, W is a pro- 
jector onto the eigenmode corresponding to the neutral 
stability with respect to spatial translations and is calcu- 
lated as the cigcnfunction £+W ~ — iwW of the adjoint 
linearized operator = D'^V^ + code + (9uf (U; Po))"^, 
sec for details [14-16]. 

For an inhomogeneity h = 9pf (U(/9, 0); po) pi(r) uni- 
form inside a disk of radius Ri, Pi(r) = H{Ri — \r\) e, 

1, where H{) is the Heaviside step function and 
e e M™, ||e|| = 1, equation (3) gives 



(4) 



Here F is the "drift force", defined as the drift velocity 
per unit inhomogeneity strength (3. For small Ri, the 
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FIG. 1: (color online) (a) Response function for a spiral wave. The spiral is visualized with a gray-scale plot of the u field. The 
yellow point indicates the spiral tip and the dashed white line shows its trajectory as the spiral rotates. Superimposed is the 
response function in term of the u component \Wu\ (red) and the v component \Wv\ (blue), (b) Enhanced visualization of one 
component of the response function. Iie{Wv) is plotted with medium gray zero (periphery), light gray positive and dark gray 
negative, (c) Drift force. The radial Fr (solid red) and azimuthal Fa (dashed blue) components of the drift force, as functions 
of the distance d, calculated by (5). The vertical dash-dotted lines show zeros of the radial component, corresponding to the 
radii at which stationary orbital movement of the spiral wave is possible. 



expression for F simplifies to 
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Fid) = / e-'^ [W{d,e)]+ apf(U(d,0);po)e 



2^' 



(5) 



We calculated the spiral wave solution U and the re- 
sponse function W using the method described in [17, 18] 
for the Barkley [19] kinetics u = {u,v), p = (a, 6, e), 
f = ifujv)'^, fu = e~'^u{l~u)(u-{v + b)/a), = u-v, 

ao = U.Y, bo = 0.1, eo = 0.02, and D ^ ^ 





This model is "excitable" , that is, it has a unique spa- 
tially uniform steady state, stable with respect to small 
perturbations [20]. At the chosen parameter values, 
the spiral wave solutions are stable [21]. Fig. l(a,b) 
shows Wu and Wy components and their location rel- 
ative to the spiral. Fig. 1(c) shows graphs of the ra- 
dial, Fr{d) = Re{F{d)) (positive for attraction) and 
azimuthal, Fa{d) ~ lm{F{d)) components of the drift 
force, for the localized inhomogeneity in parameter 6, i.e. 
e = (0,l,0)T. 

The essence of our new finding is that there is the 
change of sign of radial force Fr{d) at d = di « 3.95. 
This follows from the sign changes of W components as 
seen in fig. 1(b). For positive /3, this means attraction to 
inhomogeneity at small distances and repulsion at larger 
distances. For negative /?, however, there will be a re- 
pulsion from the inhomogeneity at o? < di and attraction 
at c? > c?i, so that d = di is a stable distance. The lat- 
ter corresponds to the drift along an orbit of radius di 
with the speed \/3Fa{di)\. There is a further root of F,. at 
d = d2 ~ 8.38; however, the corresponding value of Fa is 
very small, ^ 10~^, so no drift is easily observable there. 



Fig. 2 shows confirmation of the theoretical prediction 
of the orbital movement by direct numerical simulations 
(DNS) [18]. Panel (a) illustrates the relationship be- 
tween the DNS spiral wave solution, its tip and its in- 
stantaneous rotation centre, and panels (b) and (c) show 
the centre trajectories predicted by the theory and cal- 
culated by DNS, for Ri — 0.56 and different values of 
Sb = /3/(7ri?f). Trajectories show circular orbits, attract- 
ing for i5h < and repelling for Sb > 0, with the radius 
indistinguishable from di. 

Panels (b) and (c) illustrate two key features of or- 
bital drift: the orbiting speed depends on the strength of 
the inhomogeneity, while the radius of the orbit does not 
~ it depends only on the properties of the unperturbed 
medium. These features follow from the theory and are 
confirmed by DNS: trajectories in panel (c) have the same 
shape as in panel (b), only the spirals drift along those 
trajectories faster and in the opposite direction. 

Fig. 3 provides quantitative comparison between the- 
ory and DNS. The theoretical value of the angular ve- 
locity of the orbital movement is 17 = \Si,TTR'^Fa{di)/di\, 
implying that 17 should vary linearly with \6b\ with slope 
\nRfFa{di)/di\. Linearity of n{Sb) is indeed found in 
the DNS, remarkably up until l<5h/6o| = 0.9. The ratio 
ri/(5h in simulations is slightly smaller than the theoreti- 
cal value, due to dicretization and approximations used. 

Fig. 4 compares theory and DNS for the perturbation 
in the parameter e rather than b. There are now three 
roots of Fr{d), namely di ~ 1.97, ^2 « 3.78 and d^ « 
6.45, two of which are in the experimentally observable 
range. Roots di and d2 have alternative stability: the di- 
orbit is unstable for positive (5^ and stable for negative S^, 
and c?2-orbit is the other way round. Thus a trajectory 
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FIG. 2: (color online) Orbital movement of a spiral, (a) Simulation of a spiral wave in the presence of an inhomogeneity. The 
spiral starts near the inhomogeneity {Sb = —0.02, green disk), is repelled from it, and launches into a stable circular clockwise 
orbit [18]. The u field (red) and v field (blue) of the final spiral are shown. The preceding tip trajectory along the stable 
orbit is shown by a white line. The preceding centers rotation R are indicated both for this orbit (yellow) and the preceding 
evolution away from the inhomogeneity (blue), (b) Theoretical vector field (black arrows, nonlinearly scaled for visualization) 
and predicted trajectories (blue open circles) for the center of a spiral wave near an inhomogeneity (green disk). Actual 
trajectories for the spiral center from a DNS (red filled circles), with 5b — —0.001. Black dash-dotted circles indicate stationary 
orbits as predicted by theory. Only every 20th position of the center is shown on both theoretical and DNS trajectories, (c) 
Same for 5b = 0.003. 
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FIG. 3: (color online) Angular speed Q of orbital movement of 
the spiral wave as a function of the amplitude of the paramet- 
ric inhomogeneity Sb'. theoretical prediction vs measurements 
from direct numerical simulations. 



starting in between will enter into an orbital motion in 
any case: to the outer orbit for < and to the inner 
orbit for > 0, see fig. 4(b) and (c). 

In fig. 4(a), root di is very close to a root of Fa{d), and 
Fa{di) is very small. Thus the inner orbit is attracting for 
i5e > but the movement along it is very slow (fig. 4(c)) 
and with short observation time, it may look as if the 
spiral is attracted to any point at distance di from the 
inhomogeneity and stands still there. 

To conclude, we have reported a new type of interac- 
tion of spiral waves with localized inhomogeneities: or- 
bital movement. This new interaction has key features 
that should be observable in experiments, namely the 
orbiting speed depends on the strength of the inhomo- 
geneity, while the stationary orbit radii form a discrete 
set depending only on the properties of the unperturbed 



medium. The phenomenon is rather generic: we have 
found it in a number of other models and for more gen- 
eral shape of inhomogeneity [16]. The possibility of or- 
bital drift, related to a change of sign of an equivalent 
of Fr{d) was discussed at a speculative level in [22]; how 
often this phenomenon may occur in reality is a more 
complicated question. The equivalent of response func- 
tions calculated in [23] has a structure which suggests 
that for large-core spirals there is an infinite set of or- 
bits. In practice, orbital motion can only be observed for 
lower orbits where the orbiting speed is noticeable. The 
orbits have alternating stability, depending on the sign 
of the inhomogeneity. From this viewpoint, "pinning" as 
considered in [24] in the same model as here, appears as 
a degenerate case of orbital motion, with a zero radius. 

In certain circumstances, while an orbiting spiral may 
have the same macroscopic signature as a meandering 
one [30], the microscopic details leading to this motion are 
different. Meander, in the proper sense, is due to internal 
instabilities of a spiral wave, whereas orbital motion is 
due to inhomogeneity. E.g. in orbiting, the "meandering 
pattern" determined hy fl/uj will change depending on 
the inhomogeneity strength. 

In heart muscle, pinning of re-entrant waves of exci- 
tation has been identified as a mechanism of conversion 
of ventricular fibrillation into ventricular tachycardia [25] 
and shown to interfere with antitachycardia pacing [26]. 
Orbital movement is fundamentally more general than 
classical pinning: in pinning, the spiral is attracted at all 
distances to the inhomogeneity of a certain sign, while or- 
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FIG. 4: (color online) Orbital movement of the spiral due to perturbation in parameter e. (a) Drift force components as functions 
of distance for this inhomogeneity. The notation is the same as in fig. 1(c). (b,c) Comparison of theoretical predictions and 
simulations, for (b) = —0.001 and (c) (5e = 0.001. The notation is the same as in fig. 2(b,c). Shown are pieces of trajectories 
of the same temporal length t = . . . 500. 



bital movement occurs when repulsion at small distances 
changes to attraction at larger distances, and this can 
happen at various distances and at either sign of inho- 
mogeneity. The important practical conscquense of the 
orbital movement is that a spiral may be bound to inho- 
mogeneities of either sign, even if it is repelled from the 
inhomogeneity at small distances. 

VNB is grateful to V.I. Krinsky for inspiring discus- 
sions regarding the problem of pinning. This study has 
been supported in part by EPSRC grants EP/D074789/1 
and EP/D074746/1. 
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Calculations of the spiral wave solution U and the re- 
sponse function W. Wc used a polar grid in a disk of 
radius 12.8, discretized with 320 intervals in the radial 
direction and 128 intervals in the angular direction. 

Direct numerical simulations. Wc used forward Eu- 
ler timestepping with At = 0.00128, five-point approx- 
imation of the Laplacian with Ax = 0.08 and no-flux 
boundary conditions in a rectangular domain. The size 
of the domain was chosen big enough so that further in- 
crease did not change the behaviour. Typical sizes were 
24 X 24 and 28 x 24. Typical initial conditions were 
u{x, y, 0) = H{x — x^), v{x, y, 0) = —b+aH{y — y^:) where 
X* , were chosen depending on the desired position of 



the spiral wave. 

The tip was defined as the point where m = 0.5 and 
V = 0.25 at the given time, using bilinear interpolation 
between the grid nodes. The angle of Vu at the tip with 
respect to x axis, calculated using the same interpola- 
tion, was taken as its orientation. Positions of the centers 
were calculated by averaging the tip position during the 
time intervals when the orientation made the full circle 
(-7r,7r]. 

The angular velocity of the orbital drift fl was calcu- 
lated as 2tt/P, where P was the simulation time taken 
for one whole turn of the orbital drift, calculated to the 
nearest spiral rotation period. 



